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We present the first lattice QCD calculation of the decay constants f B and /s s with physical 
light quark masses. We use configurations generated by the MILC collaboration including the effect 
of u, d, s and c HISQ sea quarks at three lattice spacings and with three u/d quark mass values 
going down to the physical value. We use improved NRQCD for the valence b quarks. Our results 
are f B = 0.186(4) GeV, f B „ = 0.224(5) GeV, f B Jf B = 1.205(7) and M B , - Mb = 85(2) MeV, 
superseding earlier results with NRQCD b quarks. We discuss the implications of our results for the 
Standard Model rates for B( B ) — > n + /-i~ and B — > rv. 



I. INTRODUCTION 

The B and B s decay constants are key hadronic param- 
eters in the Standard Model (SM) rate for B^ — > 
and B/B s oscillations, with the B meson decay constant 
also determining the rate for B — > tv. The combination 
of experiment and theory for these processes provides im- 
portant constraints on CKM unitarity jlj and the search 
for new physics, but the strength of the constraints is 
typically limited by the errors on the hadronic parame- 
ters. 

The decay constants can only be determined accurately 
from lattice QCD calculations. Several methods have 
been developed for this [2], with errors decreasing over 
the years as calculations have improved. Here we pro- 
vide a step change in this process, giving the first results 
for /b and /s a that include physical u/d quark masses, 
obviating the need for a chiral extrapolation. As a result 
of this and other improvements described below, we have 
signficantly improved accuracy on fs a /fB over previous 
calculations. The implications of our result are discussed 
in the Conclusions. 



II. LATTICE CALCULATION 

We use eight ensembles of 'second-generation' gluon 
field configurations recently generated by the MILC col- 
laboration [U [5] , with Nf = 2 + 1 + 1 Highly Improved 
Staggered Quarks (HISQ) [B] in the sea. To control dis- 
cretisation effects, we use three lattice spacings ranging 
from 0.15 fm to 0.09 fm and light to strange mass ratios 
of mi/m s ~ 0.2,0.1,0.037. Details of the ensembles are 
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TABLE I: Details of the gauge ensembles used in this cal- 
culation. /3 is the gauge coupling, ax is the lattice spacing 
as determined by the T(2S — 15) splitting in [3], where the 
three errors are statistics, NRQCD systematics and experi- 
ment. ami,am s and am c are the sea quark masses, L x T 
gives the spatial and temporal extent of the lattices and n c f g 
is the number of configurations in each ensemble. The ensem- 
bles 1,2 and 3 will be referred to as "very coarse", 4,5 and 6 
as "coarse" and 7,8 as "fine". 

Set /3 a~c (fm) ami am s am c L x T n c f g 

1 5.8 0.1474(5)(14)(2) 0.013 0.065 0.838 16x48 1020 

2 5.8 0.1463(3)(14)(2) 0.0064 0.064 0.828 24x48 1000 

3 5.8 0.1450(3)(14)(2) 0.00235 0.0647 0.831 32x48 1000 

4 6.0 0.1219(2)(9)(2) 0.0102 0.0509 0.635 24x64 1052 

5 6.0 0.1195(3)(9)(2) 0.00507 0.0507 0.628 32x64 1000 

6 6.0 0.1189(2)(9)(2) 0.00184 0.0507 0.628 48x64 1000 

7 6.3 0.0884(3)(5)(1) 0.0074 0.037 0.440 32x96 1008 

8 6.3 0.0873(2)(5)(1) 0.0012 0.0363 0.432 64x96 621 



shown in table [Tj The lattice spacings of five of the en- 
sembles were determined using the T(2S — IS") splitting 
in [3] where details, including a discussion of the sys- 
tematic errors, can be found. The lattice spacing values 
of the additional ensembles (sets 3, 6 and 8) are deter- 
mined in the same way. The valence part of the calcu- 
lation uses lattice NonRelativistic QCD (NRQCD) [7HS] 
for the b quarks; the action is described in detail in [3J. 
It includes a number of improvements over earlier calcu- 
lations, in particular one-loop radiative corrections (be- 
yond tadpole-improvement) to most of the coefficients of 
the 0(i> 4 ) relativistic correction terms. This action has 
been shown to give excellent agreement with experiment 
in recent calculations of the bottomonium [3j [10] and B- 
meson spectrum [TT]. We are now building on previous 
calculations with the tree level NRQCD action [T2HI1] 
to extend this to B-meson decay constants. The b quark 
mass is tuned, giving the values in Table [Hi by fixing the 
spin-averaged kinetic mass with the T/rjb masses. 
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TABLE II: Parameters used for the valence quarks, arrib is 
the bare 6 quark mass in lattice units, uql is the Landau link 
value used for tadpole-improvement, and am"* 1 , am™ 1 are the 
HISQ light and strange quark masses. 



Set 


anib 


UOL 


amT 1 


amT 1 


1 


3.297 


0.8195 


0.013 


0.0641 


2 


3.263 


0.82015 


0.0064 


0.0636 


3 


3.25 


0.819467 


0.00235 


0.0628 


4 


2.66 


0.834 


0.01044 


0.0522 


5 


2.62 


0.8349 


0.00507 


0.0505 


6 


2.62 


0.834083 


0.00184 


0.0507 


7 


1.91 


0.8525 


0.0074 


0.0364 


8 


1.89 


0.851805 


0.0012 


0.0360 



TABLE III: Coefficients for the perturbative matching of the 
axial vector current (Eq. [4|. zo = po — Cio, Zi = pi — Zo, 
Z2 — pi from [15] . 



Set 


zo 


Zl 


Z2 


1 


0.024(2) 


0.024(3) 


-1.108(4) 


2 


0.022(2) 


0.024(3) 


-1.083(4) 


3 


0.022(1) 


0.024(2) 


-1.074(4) 


4 


0.006(2) 


0.007(3) 


-0.698(4) 


5 


0.001(2) 


0.007(3) 


-0.690(4) 


6 


0.001(2) 


0.007(2) 


-0.690(4) 


7 


-0.007(2) 


-0.031(4) 


-0.325(4) 


8 


-0.007(2) 


-0.031(4) 


-0.318(4) 



The HISQ valence light quark masses are taken to be 
equal to the sea mass except on set 4 where there is a 
slight discrepancy. The s quark is tuned using the rj s 
meson (M Vs = 0.6893(12) GeV [3]). Values very close to 
the sea s masses are found, meaning that partial quench- 
ing effects will be small. 

To improve the statistical precision of the correlators, 
we take £7(1) random noise sources for the valence quarks 
using the methods developed in [T3J. Along with the 
point source required for the matrix element, we include 
gaussian smearing functions for the b quark source with 
two different widths. We include 16 time sources with b 
quarks propagating both forward and backward in time 
on each configuration. We checked the statistical inde- 
pendence of results using a blocked autocorrelation func- 
tion [3J. Even on the finer physical point ensembles, the 
correlations are very small between adjacent configura- 
tions and the integrated autocorrelation time is consis- 
tent with one. 

The decay constant is defined from (0|^4o|^)qcd = 
MB q fB q , but the quantity that we extract directly from 
the amplitude of our correlator fits is &B q = yjMs q fs q ', 
we convert to f Bq at the end. For NRQCD, the full QCD 
matrix element is constructed from effective theory cur- 
rents arranged in powers of \jm\y. For Aq we consider 
the following currents, made from heavy quark \1/q and 
light quark fields 



(0) 



Jt 



J (1) - 
j — 



t(2) ^ 
J — 



2771ft 

2m& 



(1) 
(2) 
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These currents are related to the full QCD current 
through O (a s ,a s AQ C i)/m b ) by 



(A Q ) = (l + a s z )i K (4 0) ) 

+ (l + a s z 1 )(4 1) ) + a s z 2 (4 2) )^ (4) 

One-loop coefficients were calculated in [15] . Here we 
re-order the perturbation series to make the process of 



TABLE IV: Raw lattice amplitudes for B 3 and B from each 
ensemble, errors are from statistics/fitting only. a 3 ^ 2 $q and 
a 3/2^(i) arg leading amplitude and 1/mj, correction. 



Set 


a 3 / 2 3>i 0) 




a 3/2 $ (0) 


a 3/2 $ (l) 


1 


0.3720(10) 


-0.0300(3) 


0.3220(19) 


-0.0260(3) 


2 


0.3644(6) 


-0.0291(3) 


0.3093(11) 


-0.0257(8) 


3 


0.3621(16) 


-0.0288(2) 


0.2986(17) 


-0.0237(4) 


4 


0.2733(4) 


-0.0234(2) 


0.2373(9) 


-0.0197(4) 


5 


0.2679(3) 


-0.0234(1) 


0.2272(7) 


-0.0197(3) 


6 


0.2653(2) 


-0.0229(1) 


0.2193(8) 


-0.0194(3) 


7 


0.1747(3) 


-0.0170(1) 


0.1525(8) 


-0.0146(6) 


8 


0.1694(3) 


-0.0167(0) 


0.1386(5) 


-0.0136(1) 



renormalisation clearer. The Zi depend on am& and are 



given in Table III for the range of masses needed here. 
We see that the one-loop renormalisation of the tree-level 
current, Jq ^ + Jq , is tiny [3_T]. zq includes the effect of 
(o) „„A tW 



mixing between Jq and J{ at one-loop. We evaluate 
the renormalisation of Eq. [4] using a s in the V-scheme 
at scale q = 2/ a. Values for a s are obtained by running 
down from a!p(M z ) = 0.1184 [16] and range from 0.285 
to 0.314. 



III. RESULTS 



(o) 



We fit heavy-light meson correlators with both Jq 

and Jg 1 " 1 operators at the sink simultaneously using a 
multi-exponential Bayesian fitting procedure [17 j. The 



TABLE V: Raw lattice energies from each ensemble, errors 
are from statistics/fitting only. aM^ are the pion masses used 
in the chiral fits, aE(B s ) and aE(B) are the energies of the 
B s and B meson. Results on sets 3, 6 and 8 are new, others 
are given in |llj . 



Set 




aE(B s ) 


aE(B) 


3 


0.10171(4) 


0.6067(7) 


0.5439(12) 


6 


0.08154(2) 


0.5158(1) 


0.4649(6) 


8 


0.05718(1) 


0.4025(2) 


0.3638(5) 
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B and B s are fit separately; priors used in the fit are 
described in The amp litudes and energies from the 
fits are given in Tables 



IV 



andM a 3 / 2 <I>g ' is the matrix 

element of the leading current J^f 1 and a 3 / 2 ^ 1 ' that of 

Jq 1 ' and Jq 2 \ whose matrix elements are equal at zero 
meson momentum. Notice that the statistical errors in 
$ do not increase on the physical point lattices, because 
they have such large volumes. 

We take two approaches to the analysis. The first is 
to perform a simultaneous chiral fit to all our results for 

$ s , <f> s /<I> and M Bs - M B using SU{2) chiral pertur- 
bation theory. The second is to study only the physical 
u/d mass results as a function of lattice spacing. 

For the chiral analysis we use the same formula and 
priors for Mg B — M B as in [11]. Pion masses used in 
the fits are listed in Table [V] and the chiral logarithms, 
l(M%), include the finite volume corrections computed 
in [TH] which have negligible effect on the fit. For the 
decay constants the chiral formulas, including analytic 
terms up to M 2 and the leading logarithmic behaviour, 
are (see e.g. [!§]): 



$ s0 (i.o + M#/A 

$ ( 1.0 



X' 

Ml 1 



3 5 2 



2A 2 



AMD 



(5) 
(6) 



The coefficients of the analytic terms b s , bi are given 
priors 0.0(1.0) and $o,$sO have 0.5(5). To allow for 
discretisation errors each fit formula is multiplied by 
(1.0 + di(Aa) 2 + d 2 (Aa) 4 ), with A = 0.4 GeV. We ex- 
pect discretisation effects to be very similar for $ and $ s 
and so we take the di to be the same, but differing from 
the di used in the M Bs — M B fit. Since all actions used 
here are accurate through a 2 at tree-level, the prior on 
d\ is taken to be 0.0(3) whereas c?2 is 0.0(1.0). The di are 
allowed to have mild m& dependence as in [TT] . The ratio 
$ s /$ is allowed additional light quark mass dependent 
discretisation errors that could arise, for example, from 
staggered taste-splittings. 



Error % 



3>B a /3>£ 



Me 



M E 



<1>i 



4>l 



EM: 


0.0 


1.2 


0.0 


0.0 


a dependence: 


0.01 


0.9 


0.7 


0.7 


chiral: 


0.01 


0.2 


0.05 


0.05 


g: 


0.01 


0.1 


0.0 


0.0 


stat /scale: 


0.30 


1.2 


1.1 


1.1 


operator: 


0.0 


0.0 


1.4 


1.4 


relativistic: 


0.5 


0.5 


1.0 


1.0 


total: 


0.6 


2.0 


2.0 


2.1 



TABLE VI: Full error budget from the chiral fit as a per- 
centage of the final answer. 

The results of the decay constant chiral fits are plot- 
ted in Figs, [I] and [2] Extrapolating to the physical 
point appropriate to m; = (m u + m^)/2 in the absence 
of electromagnetism, i.e. M n — M w o, we find $s s = 
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FIG. 1: Fit to the decay constant ratio <&b s /$b- The fit 
result is shown in grey and errors include statistics, and chi- 
ral/continuum fitting. 
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FIG. 2: Fit to the decay constants $s s and $_g. Errors on the 
data points include statistics/scale only. The fit error, in grey, 
includes chiral/continuum fitting and perturbative errors. 



0.520(11) GeV 3 / 2 , $ B = 0.428(9) GeV 3 / 2 , $ Bs /*b = 
1.215(7). For M Bs - M B we obtain 86(1) MeV, in agree- 
ment with the result of |11) . 

Figs [3] and [4] show the results of fitting M Bs — M B 
and decay constants from the physical point ensembles 
only, and allowing only the mass dependent discretisation 
terms above. The results are <&_b s = 0.515(8) GeV 3 / 2 , 
$ B = 0.424(7) GeV 3 / 2 , $s 3 /$s = 1.216(7) and M Bs - 
M B = 87(1) MeV. Results and errors agree well between 
the two methods and we take the central values from the 
chiral fit as this allows us to interpolate to the correct 
pion mass. 

Our error budget is given in Table |VI| The errors that 
are estimated directly from the chiral/continuum fit are 
those from statistics, the lattice spacing and g and other 
chiral fit parameters. The two remaining sources of error 
in the decay constant are missing higher order corrections 
in the operator matching and relativistic corrections to 
the current. We estimate the operator matching error by 
allowing in our fits for an am&-dependent a 2 correction to 
the renormalisation in Eq. [4] with prior on the coefficient 
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FIG. 3: Fit to the mass difference Mb b — Mb on the three 
physical point ensembles only. Errors on data points include 
statistics and scale, the fit error is shown in grey. An elec- 
tromagnetic correction of -1(1) MeV has been applied to the 
lattice results and the fit to allow comparison with experi- 
ment. 
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FIG. 4: Fit to the decay constants <&b b and $s on the 
three physical point ensembles only. Errors on the data 
points include statistics/scale only. The fit error includes chi- 
ral/continuum fitting and perturbative errors. 



of 0.0(2) i.e. ten times the size of the one-loop correction, 
zo- This error cancels in the ratio f Bs /f B - We also al- 
low for a 2 corrections multiplying Jq 1 ' 2 ^ with coefficient 

0.0(1.0). The matrix element of is about 10% of J^ 0) 
from Table IV Missing current corrections at the next or- 
der in l/m b will be of size (Aqcd/w&) 2 — 0.01 which we 
take as an error. Finally, we estimated in [TT] that to 
correct for missing electromagnetic effects, M Bs — M B 
should be shifted by -1(1) MeV. 

Using the PDG masses M B , = {M B o + M B ±)/2 = 
5.27942(12) GeV and M Bs = 5.36668(24) GeV [5D] to 
convert $„ to Jb our final results are: 



fs 

!b, 

/b, / Ib 
M Bs - M B 



0.186(4) GeV 
0.224(5) GeV 
1.205(7) 
85(2) MeV. 



(7) 



For the B meson decay constant we need to distinguish 
between fs d and f Bu . Since sea quark mass effects are 
much smaller than valence mass effects we simply do 
this by extrapolating $b, and $s to values of M 2 cor- 
responding to fictitious mesons made purely of u or d 
quarks using m u /md = 0.48(10) [20]. This gives: 



W/b+ - 1-217(8) 
f B + = 0.184(4) GeV 



IbJIbo = 1.194(7) 

f B o = 0.188(4) GeV (8) 



IV. CONCLUSIONS 

Our results agree with but improve substantially on 
two earlier results using nonrelativistic approaches for the 
b quark and multiple lattice spacing values on Nf = 2 + 1 
ensembles using asqtad sea quarks. These were: f B , — 
228(10) MeV, f B Jf B = 1.188(18) (NRQCD/HISQ) JH] 
and f Bs = 242.0(9.5) MeV and f B Jf B + = 1.229(26) 
(Fermilab/asqtad) [2T]. We also agree well (within the 
2% errors) with a previous result for f Bs of 225(4) MeV 
obtained using a relativistic (HISQ) approach to b quarks 
on very fine Nf = 2 + 1 lattices [22 . Our simultaneous 
determination of M Bs — M B to 2% agrees with experi- 
ment (87.4(3) MeV 

We can determine new lattice 'world-average' error- 
weighted values by combining our results in Eq. [7] with 
the independent results of [2T] and [25] since effects from 
c sea quarks, which they do not include, should be neg- 
ligible [52] . The world averages are then: f Bs — 225(3) 
MeV and f B Jf B + = 1.218(8) giving f B + = 185(3) MeV. 

These allow for significant improvements in predictions 
for SM rates. For example, updating [24 with the world- 
average for f Bg above and our result for f B o (Eq. [8]) wc 
obtain: 

Br(B s -> /i + /i _ ) = 3.17+0.15 + 0.09 x 10~ 9 
Br{B d -> fi+fj,-) = 1.05 + 0.05 + 0.05 x 10^ 10 (9) 

where the second error from f B has been halved and is 
no longer larger than other sources of error such as V t * b V tq . 
Note that this is the flavor-averaged branching fraction 
at t = 0; the time- integrated result would be increased 
by 10% in the B s case (to 3.47(19) x 10" 9 ) to allow for 
the width difference of the two eigenstates [25] [2j)j . The 
current experimental results |27) for B s — > agree 
with this prediction. 

From the world-average f B + above we also obtain the 
Standard Model rate: 



1 



IK, 



r Br(B + ->■ tv) = 6.05(20), 



(10) 



with 3% accuracy. Calculations of matrix elements for 
B s /B mixing with physical u/d quarks are now under- 
way. 
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